% Atal, Fang, Karlsson, and Ziebarth
% this code generates main results in the paper

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    MAIN CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
addpath functions
close all

% parameters:

T = 70;
% discount rate (from institutional details)
rho = 0.966;
% interest rate:
r = 1/rho-1;
% CARA parameter (GHHW)
gamma = 0.0004;
% number of individuals for simulation
nInd = 500000;

% allow for mortality 1/yes 0/no -- should always be set to 1
mortality = 1;
   
% set seed
rng(1234567,'twister')

% replicate = 1 : uses draws used in paper
% replicate = 0 : generate new draws
replicate = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% CHOOSE FILES AND OUTPUT FOLDER  %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Income file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    income_file='../../data/processed/IncomeProcess/trends_v2.csv';

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % E(m|lambda,age) 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    data = importdata('../../data/processed/HealthProcess/mcere7A.csv');    
    
    val=data.data;
    nms = data.textdata;
    for ci=1:size(val,2)
    eval( [ nms{ci}, ' = val(:,ci);' ] ); % name the column
    end
    
    % replicate so that the format is consistent with the versions
    % that allow for the EMu to differ across groups

    mu= repmat(mu,11,1);
    agegr = repmat([0:10]',7,1);
    state = repmat(state,11,1);
    
    agegr= agegr+1;
    agegr = sort(agegr);
    emu_data = dataset(agegr,mu,state);
    emu_data = set(emu_data,'VarNames',{'agegr','mu','state'});
    
    %put in matrix form
    
    MU = [];
    for t = 1:(T)
        
        if (t>=51)
          g = 11;
        else
          g = floor((t-1)/5)+1;
        end

        % note this is transition conditional on not dying
        % in thousand dollars
        mut = emu_data.mu(emu_data.agegr==g);

        % add dead state
        
        mut = cat(1,mut,0);
        MU = cat(3,MU,mut);
    end
    
    % put in thousand dollars
    
    MU = MU/1000;
    
    % number of lambda partitions
    % without considering dead state
    age = [25:(25+T-1)]';
    Nstates = size(MU,1)-1;

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    % transition probabilities
    %%%%%%%%%%%%%%%%%%%%%%%%%%
   
    % transitions that depend on age
     data = importdata('../../data/processed/HealthProcess/prcere7bmnl.csv');  
     val=data.data;
              
        % every 5 years there is a group change
        % except after 75 years old (t = age+24)
        % so matrix is constant after t>=51 or age>=75
        % note this is transition conditional on not dying - so dividing by prob of dying
        
     CAPH = val(val(:,1)==0,3:end);
     d = cat(2,zeros(1,Nstates),1);
     CAPH = cat(1,CAPH,d);
        
            if mortality ==0
                CAPH = CAPH./(1-repmat(CAPH(:,end),1,Nstates+1));
                CAPH(:,end) = zeros(Nstates+1,1);
                CAPH(end,:) = [zeros(1,Nstates) 1];
            end
        
        
            for t = 2:(T-1)
                if (t>=51)
                    g = 11;
                else
                g = floor((t-1)/5)+1;
                end

                % note this is transition conditional on not dying
                CAPHt = val(val(:,1)==(g-1),3:end);
                CAPHt = cat(1,CAPHt,d);
                
                    if mortality ==0
                    CAPHt = CAPHt./(1-repmat(CAPHt(:,end),1,Nstates+1));
                    CAPHt(:,end) = zeros(Nstates+1,1);
                    CAPHt(end,:) = [zeros(1,Nstates) 1];
                    end
               
                CAPH = cat(3,CAPH,CAPHt);
            end
            
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                        
    %starting probability: capHtmod.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                        
    capHt = CAPH(1,:,1);
    capHtmod = capHt;
    
    % eliminate the possibility of starting dead, put all in 7; and then 
        % eliminate the possiblity of 7 and distribute equally
        capHtmod_c = capHtmod;
        capHtmod_c(1,8)=0;
        capHtmod_c(1,:)=capHtmod_c(1,:)/sum(capHtmod_c(1,:));
        capHtmod_c(1,7)=0;
        capHtmod_c(1,:)=capHtmod_c(1,:)/sum(capHtmod_c(1,:));
    
   % each scenario
        capHtmod_eye = eye(8);
        capHtmod_1 = capHtmod;
        capHtmod_1(1,:)=capHtmod_eye(1,:);
        capHtmod_2 = capHtmod;
        capHtmod_2(1,:)=capHtmod_eye(2,:);
        capHtmod_3 = capHtmod;
        capHtmod_3(1,:)=capHtmod_eye(3,:);
        capHtmod_4 = capHtmod;
        capHtmod_4(1,:)=capHtmod_eye(4,:);
        capHtmod_5 = capHtmod;
        capHtmod_5(1,:)=capHtmod_eye(5,:);        
        capHtmod_6 = capHtmod;
        capHtmod_6(1,:)=capHtmod_eye(6,:);    
        capHtmod_7 = capHtmod;
        capHtmod_7(1,:)=capHtmod_eye(7,:);  
        
  % save main starting probs in the paper
  save('./results/capHtmod_c.mat','capHtmod_c');
  % save expected claims
  save('./results/MU.mat','MU');

  for i = 1:7
      filecapH = sprintf('./results/capHtmod_%d.mat',i);
      capHtmod = sprintf('capHtmod_%d',i);
      save(filecapH,capHtmod);
  end    
  
%%%%%%%%%%%%%%%%%%%%%
% main simulations  
%%%%%%%%%%%%%%%%%%%%%

% fix model parameters

model_parameters.nInd = nInd;
model_parameters.CAPH = CAPH;
model_parameters.MU  = MU;
model_parameters.rho = rho;
model_parameters.r = r;
model_parameters.income_file = income_file;
model_parameters.gamma = gamma;
save('./results/model_parameters_main.mat','model_parameters');
 
load('./results/model_parameters_main.mat')
load('./results/capHtmod_c.mat')
model_parameters.outputfile ='./results/allstates_main.mat';
model_parameters.capHtmod= capHtmod_c;

if replicate ==1
    % use draws in paper
    load('./results/model_parameters_main.mat')
    load('./results/capHtmod_c.mat')
    load('./results/iX_submission/iX_main.mat')
    model_parameters.outputfile ='./results/allstates_main.mat';
    model_parameters.capHtmod= capHtmod_c;
    model_output_h = solve_main(model_parameters,iX);
    for i=1:7
       filecapH = sprintf('./results/capHtmod_%d.mat',i);
       fileiX = sprintf('./results/iX_submission/iX_%d.mat',i);
       model_parameters.outputfile =sprintf('./results/state%d.mat',i);
       load(filecapH);
       load(fileiX);
       model_output = solve_main(model_parameters,iX);
    end

else
    % generate new draws
    model_output_h = solve_main(model_parameters);
    fileiX = './results/iX_out/iX_main.mat';
    iX = model_output_h.iX;
    save(fileiX,'iX');
    for i=1:7
       filecapH = sprintf('./results/capHtmod_%d.mat',i);
       load(filecapH);
       capHtmod_name = ['capHtmod_',num2str(i)];
       model_parameters.capHtmod= eval(capHtmod_name);
       model_parameters.outputfile =sprintf('./results/state%d.mat',i);
       model_output = solve_main(model_parameters);
       % save draws for future replication
       iX = model_output.iX;
       fileiX = sprintf('./results/iX_out/iX_%d.mat',i);
       save(fileiX,'iX');
    end
end



